########################################################################
#   Project: The Media Smells like Sulfur!!! Leaders and Verbal        #
#            Attacks against the Fourth Estate in Unconsolidated       #
#            Democracies                                               #
#   Authors: JA Solis  & I Sagarzazu                                   #
#   Task: Replication Materials for Political Communication article    #
#   R version: 3.6.0                                                   #
#   Last updated: 8/14/2019                                            #
########################################################################

#Clear environment
rm(list=ls())

#set working directory
setwd("[insert working directory]")

#Load package for output
library("texreg")

####################################################################
# NOTE: Figures and Tables in manuscript made in RStudio and Stata #
####################################################################

#Load data#
ap <- read.csv("polcom_rep_ap.csv")

##############################################
# Figure 1: Please see .do file to replicate #
##############################################

###########################################
# Figure 2: Plot Quarterly verbal Attacks #
###########################################

#Make object of ts variable of interest
quarterly <- ts(ap$attack_sent, frequency=4, start=c(1999,2))

#Plot series, raw data
plot(quarterly,
     xlab="Year", ylab="Count of Verbal Attacks")


#################################################
# Figure 3: Please see .do file for replication #
#################################################

#########################################
# Table 3: Box-Taio Intervention Models #
#########################################

#NOTE: Identification process yielded arima(0,1,1). Code detailing identification process
#      avaialble upon request

#Fit ARIMA model and extract residuals
fit<-arima(quarterly, order= c(0,1,1))
attacks_resid<-residuals(fit)

#Box-Tiao Intervention Model: Regress residuals on theoretically selected pulses (and controls)#

#-#-#-#-#
#Model 1#
#-#-#-#-#

noninst_inst <- lm(attacks_resid ~ non_inst + inst + approval + oil, data=ap)
summary(noninst_inst)
AIC(noninst_inst)
#486.8856

#Create LaTex table:
texreg(list(noninst_inst), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")

#-#-#-#-#
#Model 2#
#-#-#-#-#

bevent <- lm(attacks_resid ~ cre1999q + elect2000q + prtst2002q + coup2002q 
             + strk2002q + ref2004q + elect2006q + cre2007q 
             + cre2009q +approval + oil , data=ap)
summary(bevent)
AIC(bevent)
#491.1601

#Create LaTex table:
texreg(list(bevent), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")


#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
#################################
###########  Appendix ###########
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
#################################

#################################
# Figure 1: ACF and PACF Graphs #
#################################

#ACF
acf(quarterly,main="")

#PACF
pacf(quarterly,main="")

###########################
# Table 1: Ljung-Box Test #
###########################

#Ljung-Box test (check for no serial autocorrelation among lags)
#Null Hypothesis=autcorrelation exists; alpha=.05

Box.test(attacks_resid, lag=5, type = "Ljung-Box")
#p=0.887

Box.test(attacks_resid, lag=10, type = "Ljung-Box")
#p=0.6569

Box.test(attacks_resid, lag=15, type = "Ljung-Box")
#p=0.7819

Box.test(attacks_resid, lag=20, type = "Ljung-Box")
#p=0.6556

Box.test(attacks_resid, lag=25, type = "Ljung-Box")
#p=0.7067

Box.test(attacks_resid, lag=30, type = "Ljung-Box")
#p=0.6816

Box.test(attacks_resid, lag=35, type = "Ljung-Box")
#p=0.8054

Box.test(attacks_resid, lag=40, type = "Ljung-Box")
#p=0.8423

############################################
# Figure 2: Raw vs. ARIMA(0,1,1) Residuals #
############################################

ts.plot(quarterly, attacks_resid, main="Qrtly Verbal Attacks vs. ARIMA (0,1,1) Filtered Residuals",
        xlab="Year", ylab="",gpars = list(col = c("black", "red"))) 

################################################
# Table 2: Please see .do file for replication #
################################################

#################################################
# Figure 3: Histogram of ARIMA(0,1,1) Residuals #
#################################################

hist(attacks_resid, xlab="ARIMA filtered residuals", main="")

#########################################################
# Table 3: Summary Statistics of ARIMA(0,1,1) Residuals #
#########################################################

summary(attacks_resid)
#   Min.   1st Qu.   Median    Mean    3rd Qu.   Max. 
#-58.9638 -20.6620  -4.7881  -0.8176  15.3048  95.6277
sd(attacks_resid)
#29.03375

#################################################
# Figure 4: Please see .do file for replication #
#################################################

###############################################
# Table 4: B-T Analysis, include 2007 Protest #
###############################################

#-#-#-#-#
#Model 1#
#-#-#-#-#

rctvm1 <- lm(attacks_resid ~ non_inst_rctv + inst_rctv + approval + oil, data=ap)
summary(rctvm1)
AIC(rctvm1)
#493.1455

#Create LaTex table:
texreg(list(rctvm1), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")

#-#-#-#-#
#Model 2#
#-#-#-#-#

rctvm2 <- lm(attacks_resid ~ cre1999q + elect2000q + prtst2002q + coup2002q 
             + strk2002q + ref2004q + elect2006q + cre2007q + prtst2007q
             + cre2009q +approval + oil , data=ap)
summary(rctvm2)
AIC(rctvm2)
#492.8599

#Create LaTex table:
texreg(list(rctvm2), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")

########################################
# Table 5: w/ GDP & Non-verbal Attacks #
########################################

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
##-#-#-#-#-#-## Pooled Events ##-#-#-#-#-#-##
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

########################
##########GDP###########
########################

gdp <- lm(attacks_resid ~ non_inst + inst + approval + oil + gdp_qtr, data=ap)
summary(gdp)
AIC(gdp)
#486.5765

#########################
#########################
##### Non-Vbl Atks-1 ####
#########################
#########################

nv_atks <- lm(attacks_resid ~ non_inst + inst + approval + oil + censor_fill + harass_fill, data=ap)
summary(nv_atks)
AIC(nv_atks)
#490.3417

#Create LaTex table:
texreg(list(gdp, nv_atks), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")

########################################
# Table 6: w/ GDP & Non-verbal Attacks #
########################################

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
##-#-#-#-#-#-##Individual Events##-#-#-#-#-#-##
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

########################
##########GDP###########
########################

gdp_ind <- lm(attacks_resid ~ cre1999q + elect2000q + prtst2002q + coup2002q 
                + strk2002q + ref2004q + elect2006q + cre2007q 
                + cre2009q  + approval + oil + gdp_qtr, data=ap)
summary(gdp_ind)
AIC(gdp_ind)
#491.861

########################
#######Atks-1###########
########################

nv_atks_ind <- lm(attacks_resid ~ cre1999q + elect2000q + prtst2002q + coup2002q 
            + strk2002q + ref2004q + elect2006q + cre2007q 
            + cre2009q + approval + oil + censor_fill + harass_fill, data=ap)
summary(nv_atks_ind)
AIC(nv_atks_ind)
#494.9151

#Create LaTex table:
texreg(list(gdp_ind, nv_atks_ind), dcolumn = TRUE, booktabs = TRUE,
       use.packages = FALSE, label = "tab:3", caption = "", include.bic = FALSE,
       float.pos = "hb!")

#################################################
# Figure 5: Please see .do file for replication #
#################################################


##########################################################
# Figure 6: Compare verbal attacks to episodes sentences #
##########################################################

#Note: Use `ap_episodes.csv' for this analysis

#Load new dataset
scatter <- read.csv("ap_episodes.csv")

#Create objects for plot
vattacks <- scatter$attack_sent
sent <- scatter$sentences

plot(sent, vattacks, main = "",
     xlab = "Sentences in Episodes", ylab =  "Verbal Attacks against Media",
     pch = 19, frame = FALSE)
abline(lm(vattacks ~ sent, data = scatter), col = "blue")


                                         #########
                                         ## END ##
                                         #########
